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Abstract 
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I. INTRODUCTION 



The generation problem of gravitational waves (GWs) for inspiralling compact binaries has 
been completed at the third post-Newtonian (3PN) order both for the equation of motion of the 
binary and for its far-zone radiation field. The computations of the 3PN accurate equations of mo- 
tion (EOM) and mass quadrupole moment were technically more involved than the corresponding 
2PN cases due to the issues related to the ambiguities of self-field regularisation using Riesz or 
Hadamard regularisations [0, 0, 0, 0, 0] . A deeper understanding of the cause of these ambigu- 
ities and the use of the efficient dimensional regularisation scheme was crucial to the resolution 
of this problem [0, 0, 0, 0]. The 3.5PN phasing of inspiralling compact binaries (ICBs) moving 
in quasi-circular orbits is now complete and available for use in GW data analysis [0,0,0]. This 
is timely since prototype binary GW sources for laser interferometer detectors are neutron star or 
black hole binaries close to their merger phase and consequently moving in quasi-circular orbits. 
However, astrophysical paradigms do exist that result in binaries with nonzero eccentricity in the 
sensitive bandwidth of both terrestrial and space-based GW detectors ifiol fll. 12, 13, HI 151. 

More precisely, there currently exists a variety of astrophysical scenarios that can produce bina- 
ries which have residual (~ 0.01), moderate (~ 0.5) or even very high (~ 0.9) eccentricities close to 
merger — contrary to garden- variety ICB's mentioned earlier. For instance, the Kozai mechanism 
is one important scenario that produces eccentric binaries and involves the interaction between a 
pair of binaries in the dense cores of globular clusters 111 ill . If the mutual inclination angle of the 
inner binary is strongly tilted with respect to the outer BH, then secular Kozai resonance II 1 Oh can 
increase the eccentricity of the inner binary to large values. Numerical investigations 1151 1 have 
shown that intermediate mass BH binaries can have eccentricity of order 0.9 when they are visible 
in the LISA band. Even for the case of stellar mass BHs, the eccentricity may be about 0.1 at 10Hz 
making them possible important sources for future ground-based GW detectors such as the Ein- 
stein Telescope (ET) which optimistically would attempt to achieve a seismic cut-off frequency 
around one Hertz. In the context of the "final parsec problem" for galaxy mergers, Ref. II 1611 
pointed out that angular momentum loss to the circumbinary gas, which can provide a mechanism 
for overcoming this problem, can produce small but non zero eccentricity via interaction of the 
binary with the gas disk. The resultant eccentricity can range from 0.01 - 0.1 one week prior to 
merger depending on the binary's mass ratio and will have observable effects on the LISA signal. 
A more recent study of the scattering of stellar mass BHs in the galactic centres I117ll has found 
that more massive BHs dominate the scattering rate close to the central supermassive BH. These 
scatterings could give rise to bound binaries which will have a high eccentricity (~ 0.9) when 
they enter the LIGO band. More importantly, due to higher harmonics present in the GW signals 
from eccentric binaries, these sources can be observed to larger distances and with larger masses 



(< 700 M Q ) than for circular orbits, depending on the eccentricity IU7U l8ll. This will have impli 



cations for sources in Advanced LIGO/VIRGO and ET detectors due to their very good proposed 
low frequency sensitivity 111 811 . For a detailed discussion about various astrophysical mechanisms 
related to the eccentric orbit binaries see the introduction of Ref. lfl9ll and appendix A of [fl8ll . 

The complete (ambiguity-free and fully determined) 3PN accurate EOM and mass quadrupole 
moment for compact binaries enable one to compute the 3PN energy and angular momentum 
fluxes for inspiralling compact binaries moving in general non-circular orbits. Recently, in two 
related papers [2(1 21], we laid out the formalism and implemented the computation of the GW 
energy flux for non-circular orbits up to 3PN order. For non-circular orbits, to determine the orbital 
phasing, and the secular evolution of the orbital elements, the GW angular momentum flux needs 
to be known in addition to the energy flux. Of course, we also need the conserved center-of-mass 
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energy and angular momentum of the orbit as deduced from the EOM. 

In this paper, we compute the angular momentum flux of inspiralling compact binaries up to 
3PN order generalising earlier work by Peters [12211 at Newtonian order, extended in Ref. Il23ll at 
1PN order, in Q at 1.5PN (tails) and in [0] at 2PN. The 3PN contributions to energy and an- 
gular momentum fluxes come not only from instantaneous terms but also (non-linear) hereditary 
contributions 11261 . 12711 . We shall find that for the angular momentum flux, the hereditary contri- 
butions comprise not only the tails, tails-of-tails and tail-square terms as for the energy flux but 
also formally an interesting memory contribution at 2.5PN order. One can then average the 3PN 
energy and angular momentum fluxes over an orbit thanks to the 3PN generalized quasi-Keplerian 
parametrization of the binary's orbital motion rf28ll . Finally, we compute the secular evolution of 
orbital elements under 3PN gravitational radiation reaction (i.e. corresponding formally to 5.5PN 
terms in the EOM). Th is g eneralises the works of Peters and Mathews n2% at 2.5PN, Blanchet 
and Schafer at 3.5PN [H and 4PN 0, [U orders, and Gopakumar and Iyer [HI at 4.5PN. 
While 11231 12411 require the 1PN accurate orbital description [32], Ref. I125ll crucially employs the 
generalised 2PN quas i-Keplerian parametrization of the binary's orbital motion in ADM coordi- 
nates as given in 11331 [34], HH]. In the present case the averaging of the instantaneous terms will 
require the full 3PN generalised quasi-Keplerian representation. However, the hereditary terms 
being relatively of higher PN orders, only require the 1PN parametrisation of the motion. 

The secular evolution of orbital elements under gravitational radiation reaction provides the 
starting point for constructing templates for eccentric orbits. One of the first works in this direction 
was Ref. 113 611 which investigated the efficiency with which circular-orbit based templates would 
be able to detect an eccentric-orbit signal. To go beyond the secular evolution in the gravitational 
wave phasing one needs to include besides the averaged contribution the oscillatory terms in the 
evolution of orbital elements. Damour, Gopakumar and Iyer B7I1 discussed an analytic method 
for dealing with this issue at the leading radiation reaction order of 2.5PN, making possible the 
construction of high accuracy templates for the GW signals from ICBs in quasi-elliptical orbits. 
This was extended to 3.5PN order in Ref. II 1911 . and the problem was revisited in a more elaborate 
way in 113 811 . including the computation of the noise- weighted overlaps for different astrophysical 
situations. Further investigations would be necessary to tackle the data analysis issues especially 
if the signals have moderate or high eccentricities. Including PN corrections to higher order in the 
evolution of orbital elements is a crucial step towards this. 

With the recent advances in numerical relativity (NR) has emerged the possibility of comparing 
the NR waveforms to the PN results and exploring the regime of validity of the PN approxima- 
tion 1139L |40L I41U42II. These comparisons have been done for different source configurations such 
as nonprecessing spinning binaries [|43ll . precessing spins H4\ and, very recently, eccentric bi- 
naries j45ll . It is obviously crucial to have a very accurate PN expression for the evolution of 
the GW phase while comparing with the high accuracy NR simulations. For the circular orbit 
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case we have 3.5PN accurate expressions in phasing H, & S, |47|] and 3PN accurate ones in 
amplitude [48, 49, 5(1 Hi 52]. Notably Ref. ll45ll presented the first comparison between NR sim- 
ulations of an eccentric binary black hole system with the corresponding current best PN results. 
The simulations relate to equal-mass, nonspinning binaries but with an eccentricity e ~ 0.1 for 
about twenty GW cycles before merger, and comparison to a currently available 2PN eccentric 
binary. The work 114511 explores the parametrisation most suited for such a comparison and stands 
to improve with the 3PN model that our present paper will now provide. 

The organization of this paper is as follows. In Sec. [TT] we start with the basic expression of 
the far-zone flux of angular momentum, employ expressions relating the radiative moments to the 
source moments and decompose the angular momentum flux into its instantaneous and hereditary 
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parts. Sec.|ni]discusses the computation of the instantaneous terms in both standard harmonic co- 
ordinates and ADM coordinates. Using the 3PN quasi-Keplerian representation, Sec.lTVlcomputes 
the orbital average of the instantaneous part of the angular momentum flux in ADM coordinates. 
Sec. |V] deals with the computation of the hereditary contributions in the averaged angular momen- 
tum flux using a Fourier domain decomposition. The evolution of the orbital elements, in ADM 
coordinates, including both instantaneous and hereditary terms, is presented in Sec.|VIJ In the final 
Sec. I VIII we provide simpler explicit expressions of various inputs up to first order in e 1 needed to 
deal with phasing in the small eccentricity limit e — » 0. The paper concludes with three appen- 
dices. Appendix lAl presents an analysis of non-linear memory leading to a DC term arising from 
the dependence over the binary's (remote) past history. Appendix |B] includes tables of numerical 
values of the various "enhancement functions" appearing in the paper to facilitate comparisons 
with numerical relativity runs and use in data analysis applications. The paper concludes with Ap- 
pendix where the important equations are also presented in the modified harmonic coordinates 
for the convenience of the user. 



H. THE FAR-ZONE ANGULAR MOMENTUM FLUX 

In this section, we start from the angular momentum flux expressed in terms of the radiative 
multipole moments, use the relations connecting those radiative moments to the source moments, 
and rewrite the flux as a sum of the instantaneous terms which are functions of the retarded time, 
and hereditary terms which depend on the dynamics of the system in its entire past. The 3PN 
accurate angular momentum flux in the source's far- zone, denoted for convenience 



GW 



(2.1) 



is expressed in terms of the mass and current type radiative multipole moments in radiative coor- 
dinates [53] as 
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In the above U L and V L (with L = ■ ■ ■ k a multi-index composed of / indices) are the symmetric- 
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trace-free (STF) mass and current type radiative multipole moments respectively, and U[ and V 
denote their p th time derivatives. The moments are functions of retarded time T R = T - R/c in the 
radiative coordinates (T, X), with R = |X| the distance of the source; e,^ is the usual Levi-Civita 
symbol such that £123 = +1; the shorthand 0(n) denotes a PN remainder of order of 0(c~ n ). 



A. Radiative moments in terms of source moments 



Using the MPM formalism II54L 15511 . the radiative moments in Eq. (12.21) can be computed in 
terms of the source moments to an accuracy sufficient for the computation of the angular momen- 
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turn flux up to 3PN. One must compute the mass type radiative quadrupole to 3PN accuracy, 
mass octupole Ujjk and current quadrupole Vjy to 2PN, mass hexadecupole Uijk m and current oc- 
tupole Vijk to 1PN, and finally f/^ mn and V/y^n to Newtonian order only. The relations connecting 
the radiative moments U L and V L to the corresponding source moments L and J L (and also to 
the so-called gauge moments W L , X L , Y L and Z L ) are now given l54|,|55|]. For mass type 

moments we have (the brackets <> denoting the STF projection) 
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The Zl's and 7 L 's are the STF mass and current-type source moments; M = I is the total mass 
monopole which is conserved; W is the monopole corresponding to one type of the gauge mo- 
ments, i.e. W L . For the current-type moments we find 
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V ijk (T R ) = J®(T R ) + OQ), 
V iM (T R ) = J^ kl (T R ) + 0(3). 
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The radiative moments have two distinct contributions: "instantaneous", which is a snapshot 
function of the retarded instant T R only; and "hereditary", which depends on the dynamics of the 
system in its entire past V < T R . The parameter To appearing in the logarithms of Eqs. (12.31) 



and (12.41) is a freely specifiable constant time scale, entering the relation between the retarded time 
T R = T-R/c in radiative coordinates and the corresponding time ? H -r H /c in harmonic coordinates 
(where r H is the distance of the source in harmonic coordinates). Posing r = cr we have 
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We choose the constant r scaling the logarithm to match with the choice made in the computation 
of tails-of-tails in Il26n . 
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B. Structure of the 3PN angular momentum flux 



One can schematically split the total contribution to the angular momentum flux as the sum of 
the instantaneous and hereditary terms, 



where the instantantaneous terms are given by [not specifying the obvious 0(n)] 
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Using (12.31) — (12.41) in Eq. (12.21) the hereditary part can be further decomposed as 
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where the quadratic tail integrals are given by 1 
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1 Ref. 112411 contains a typographical mistake at 1.5PN in Eq. (83). The first term in Eq. d2.9l l below contai ning 
"7 (2) (7r)/ (5) ( V)" is missing. This has been independently pointed out recently in 15611 . However, the results in 124 ] 
are correct and do take this term into account. Ref. 12511 which only quotes 12411 also contains this typo. 
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These tail and tail-of-tail integrals are similar to those occurring in the energy flux piQ . However 
we have also the non-linear memory integral 
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Recall that the memory contributes to the radiative quadrupole moment Uu simply by an anti- 
derivative of source moments. It therefore becomes instantaneous when we consider the time 
derivative if:) of the radiative quadrupole moment. Hence we have incorporated in the instanta- 
neous part of the angular momentum flux [Eq. ( 12.71 )1 a term coming from the time derivative of 
the memory integral. The presence of the non-linear memory contribution ( 12.1 II ) remaining in the 
angular momentum flux is to be noticed. This is in contrast with the case of the energy flux, where 
there is no memory contribution because the memory is time differentiated therein and therefore 
becomes instantaneous (see Ref. 112 ill ). 



III. INSTANTANEOUS TERMS IN THE 3PN ANGULAR MOMENTUM FLUX 

The relevant source multipole moments needed for the computation of the angular momentum 
flux up to 3PN order are the same as for the energy flux. Hence we redirect the reader to Sees. Ill 
and IV of Ref. [20] for a detailed discussion of these moments as well as for the equation of 
motion up to 3PN order (see Ref. |[57ll ) which is required when differentiating the source moments. 



Notice that the mass quadrupole moment 7, ; was available in [5] and was used in 112 111 . Using all 
these source moments, and the EOM at 3PN order, it is possible to compute the different PN 
contributions to the instantaneous part (12.71 ) of the angular momentum flux. 

The prominent application of the present computation will be discussed in Sec. [VTl where the 
evolution of the orbital elements under gravitational radiation reaction will be investigated to 3PN 
order. This will be based on the solution of the motion in the form of the quasi-Keplerian repre- 
sentation of the orbit. The latter can be written up to 3PN order in ADM coordinates (and also in 
"modified" harmonic coordinates 1I20I1 ) but not in the standard harmonic coordinate system, due to 
the presence of gauge-dependent logarithms arising at 3PN order. We shall therefore present our 
results in the (standard) harmonic coordinate system and also in ADM coordinates, in which we 
give the evolution of the orbital elements (see Sec. [VTb . 
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A. Instantaneous flux in standard harmonic coordinates 



The standard harmonic (SH) coordinate system is the one in which the original computations 



of the 3PN center-of-mass equations of motion Il57h and 3PN source quadrupole moment [5] were 



given. It is known that the standard harmonic coordinate system develops some gauge-dependent 
logarithmic terms at 3PN order. This particular coordinate system is referred to as standard in 
order to distinguish it from the modified harmonic (MH) coordinate system, where the logarithms 
at 3PN order have been gauged away. Both the MH coordinates and the ADM coordinates are 
suitable for a 3PN quasi-Keplerian parametrization of the motion; they were reviewed and used in 
the computation of the energy flux in Ref. rf20ll . 

The angular momentum flux will be orthogonal to the orbital plane, and aligned with the orbital 
angular momentum. We introduce the unit direction along the orbital angular momentum, 

(3.1) 

where x* and v k = dx k /dt are the relative separation and velocity of the two particles. Here 
<p = d(f/dt, where <p denotes the orbital phase 2 , and is linked to the orbital velocity by v 2 = r 2 + r 2 <p 2 
where r = dr/dt. Posing 

= A'&nst , (3.2) 

we look for the 3PN expansion 

£?inst = + £?1PN + ^2PN + £?2.5PN + £?3PN + 0(7) , (3.3) 

where as usual the Newtonian piece really corresponds to the dominant radiation reaction at 2.5PN 
order in the equations of motion. The results in SH coordinates then read 3 
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below, namely <p = co + n dW/d(, where W{() is periodic in I — n(t - tp) with period 2n, hence periodic in time with 



period P (see e.g. Sec. II. A in 0210 ). We can call <p the "instantaneous" angular frequency. 
3 Mass parameters are the total mass m = ni\ + mi, the reduced mass [i = mimi/m and the symmetric mass ratio 
v = fi/m which is such that < v < 1/4. 
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(3.4e) 



We reproduce the terms computed earlier 11221 |23L 12511 in the angular momentum flux up to 2PN 
order. As can be seen from above, the 3PN terms contain two kinds of logarithms, ln(r/r ) and 
ln(r/r' Q ). The logarithms ln(r/rQ) are specific to the standard harmonic coordinate system. We 
discuss later in more detail these different kinds of logarithmic terms. 



B. Instantaneous flux in ADM coordinates 

Since many related numerical relativity studies are in ADM-type coordinates, we shall present 
the applications in later sections of this paper in ADM coordinates. Going from SH to ADM coor- 
dinates removes the gauge-dependent logarithms ln(r/r' ) which are not very convenient to handle 



in numerical calculations. We shall thus re-express the instantaneous flux (I3.3I) - (I3.4I) in terms of 
the variables in ADM coordinates. This requires the use of the so-called contact transformation 
of these variables {x l , v l and r), linking the standard harmonic coordinates to the ADM ones. We 
refer to lf2Qh (see Sec. VLB) for the details of that transformation. The contact transformation 
equation will depend on some ln(r/rQ)'s and that dependence will be such that the flux in ADM 
coordinates becomes independent of those logarithms ln(r/rQ), revealing the gauge nature of the 
constant r' Q . However, as we shall comment below, the other logarithms of type ln(r/r ) will re- 
main in the instantaneous part of the flux but will be cancelled by related contributions coming 
from the tail integrals (and more precisely the tails-of-tails); indeed, see the constant To = r /c in 
Eqs. (I2~9l)- (ITT01) . 

The angular momentum flux in ADM coordinates admits the same type of PN expansion 
as (13 ,3b . Since the transformation between the SH and ADM coordinate systems starts at the 
2PN order, only the 2PN and 3PN terms in the flux get modified and we now list these two terms, 
labelled by "ADM" to remember that all variables therein correspond to ADM coordinates. Let 
us recall that the transformation formulas between the standard harmonic and ADM coordinates 



which were given in [20] concerned only the conservative part of the dynamics. We thus give here 
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only the 2PN and 3PN terms in the angular-momentum flux following these transformations, 
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(3.5b) 



C. 2.5PN terms in the fluxes in ADM coordinates 



In the previous section, we have taken into account the contact transformation involving "con- 
servative" orders up to 3PN required to go from the standard harmonic coordinates to the ADM 
coordinates. However, there still remains the possible change of gauge in the radiation reaction 
(dissipative) terms at order 2.5PN. We discuss these terms here since the corresponding transfor- 
mation law was not given in Ref. [|20ll . 4 We recall that in the standard harmonic (SH) coordinate 
system the lowest-order dissipative part of the equations of motion, i.e. the 2.5PN acceleration 
term, is given by (with boldface letters indicating ordinary three-dimensional vectors) 



8 G 2 m 2 v 



SH 

l 2.5PN- 5 ^3 



3v 2 + 



17 Gm 





9 „Gm 




rn + 


-v 2 -3 


•} 




r 





(3.6) 



4 Actually the dissipative terms at 2.5PN order in the fluxes are not very important for the present purpose because 
they will average to zero and thus will not contribute to the balance equations. The present discussion is for 
completeness and future use. 
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However, one may choose to work in alternative radiation gauges and a convenient characterisation 
at 2.5PN has been investigated earlier in rf58l [59h (see also I160ll ). Following this work the most 
general form of the relative acceleration is specified by the two-parameter family written as, 



l 2.5PN 



8 G 2 m 2 v 
5 c 5 r 3 



(A 2 .5PN^n + 52. 5 p N v), 
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3 r 

B 2 5pn = -(2 + a)v 2 - (2 - a)— + 3(1 + a)r 2 . 

r 



(3.7a) 
(3.7b) 
(3.7c) 



The general 2.5PN gauge is parametrized by the two numerical constants a and ft. The SH gauge 
in which the acceleration is given by (13.61) corresponds to setting a = -1 and /3 = 0; the ADM 
gauge is obtained by using a = 5/3 and = 3, in which case the 2.5PN acceleration becomes 116 ill 



a 



ADM 
2.5PN 



8 G l m l v 
5 c 5 r 3 



12 v 2 + 2 

r 



\5r 2 



r n + 



11 



1 Gm 



ir - — - 



3 r 



+ 8r 2 



(3.8) 



To obtain the energy and angular momentum fluxes in a general 2.5PN gauge parametrised by 
a and /? we apply upon the SH particle's worldlines the (o^/^-dependent shift 11581. 159TI 



8 G 2 m 2 v 
15 c 5 r 



(-/?rn + (3 + 3a-2/3) v) , 



(3.9) 



such that the relative position of the bodies in a general gauge is given by x gen (0 = x SH (0 + e(t). 
The relative velocity and acceleration of the bodies are then v gen (0 = v SH (0 + de/dt and a gen (0 = 
a SH (0 + d 2 e/dt 2 . However, the effect of this shift on the acceleration is more subtle. Indeed, for 
the acceleration, one must remember that it is a functional of the position and velocity coming 
from the equation of motion and consequently the acceleration must consistently be expressed in 
terms of the position and velocity in the corresponding frame. Thus a gen (?) = a SH (0 + d 2 e/dt 2 
is correct, but here our too condensed notation for the acceleration means in fact that a gen (0 = 
a gen [x gen (0,v gen (0] in the general frame and a SH (f) = a SH [x SH (f), v SH (f)] in the SH frame. Re- 
expressing the latter relation in terms of "dummy" variables x(f) and \(t) we then get the functional 
relation a gen [x(f), v(/)] = a SH [x(?X + 8 £ a(t) where the change in the acceleration is given by 



d 2 e 
6 ^=df 



da N 



(3.10) 



Here, neglecting higher-order terms, aN denotes the ordinary Newtonian acceleration which de- 
pends on x but not on v. 

Consider now the effect of the shift on the components of the quadrupole moment (the gener- 
alization to higher moments is trivial). Under this shift the mass quadrupole moment Ifj 1 in SH 
coordinates will be changed into lfj n such that 7 ; ^ en [x gen , v gen ] = 7 ; ^ H [x SH , v SH ]. Hence, introduc- 
ing dummy variables x, v we readily obtain ^ e "[x, v] = 7^ H [x, v] + 8Jij where the change can be 
expressed in terms of the Newtonian quadrupole moment / 



(N) 



as 5Jij = 



-e k dl ( j j ) /dx k (neglecting 
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higher-order PN terms). Using 7^ = mvx^x^ (which depends only on x) we get, 5 

SJij = -2mvjcV. (3.11) 

Next, we compute the successive time derivatives of this quadrupole moment, being careful that 
when reducing the result by means of the equation of motion the acceleration is modified by the 
amount (13 . 1 OP . This order reduction starts from the second time derivative. We readily obtain 
d 2 /f ; en [x, \]/dt 2 = (ag en [x, v]d/™[x, \]/dx k + •••) + d 2 6Jij/dt 2 , where we have explicitly replaced 
the acceleration by the general-frame one (since we are evaluating in the left-hand-side the time 
derivative of the general-frame moment), and where the dots indicate other terms which do not 
need any replacement of accelerations. Hence, using the modification of the acceleration as given 
by (13.101) . we can connect the right-hand-side to what we would get for the time derivatives of the 
SH moment using the SH acceleration. This yields /ff n [x, v] = /f^fx, v] + Sja where the total 

'J 'J ■* 

change is now given by 8j]j = 8 £ a k dlfj/dx k + d 2 S e Iij/dt 2 (the dots mean the time derivatives). The 
results, which proceed in the same way for the third time-derivative, read 

d 2 

Sj'ij = 2m v x <[ 6 £ a j> + — (8J U ) , (3.1 2a) 

+ ^(SJ V ). 0.12b) 



8j ij = 2m v 



3v <! 8 f a }> + x <l 



dt 2 

d6 f a j> 



dt 



Beware that the total change in the time derivative of moments 8jij is different from the time 
derivative of the change of the moments d 2 8Jij/dt 2 . 

Finally we obtain the shifts of the energy and angular momentum fluxes. Since they are given in 
first approximation by the quadrupole formulas, for instance !F gen [x, v] = ^ I f, en [x, v] I f, e "[x, v], 

we can immediately use the previous link 7 gen [x, v] = /ff[x, v] + 8 £ Ijj to conclude that the 
associated changes in the fluxes, say !F gen [x, v] = T $li [x, v] + 8 £ T, are 

8 £ r = ^- 5 'l^ ) Sj ij , (3.13a) 

Se&i = |j Siab [C * J he ~ "ifc 8ehc] ■ (3.13b) 

Evidently the factors of the gauge modified multipole moments are Newtonian at this order and 
straightforwardly computed. Explicit results for the 2.5PN terms in the energy and angular mo- 



5 Some signs should be corrected in Ref. II20I1 : namely, Eq. (6.4) should have a minus sign, Eqs. (6.5) and (6.7) have 
plus signs, and the second term of Eq. (6.6) should have a plus sign. The results in Ref. 12011 . notably Eq. (6.8), are 
unchanged. 
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mentum fluxes in the (a,/?)-dependent gauge are then 6 
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(3.14b) 



These expressions could be of use in the discussion of non-quasi-circular effects required to match 
high accuracy PN waveforms to those of numerical relativity in the effective-one-body formal- 
ism 11621 6311 . By setting a = -1 and = one recovers the SH coordinates results of Eq. (5. 2d) 
of Ref. [|20l] and (|3.4dt of this paper. On the other hand, for a = 5/3 and fi = 3 we get the ADM 
results given by 
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(3.15b) 



From now on our calculations will be in ADM coordinates, so we henceforth suppress all indica- 
tions regarding the ADM coordinate system. 



IV. ORBITAL AVERAGE OF THE INSTANTANEOUS ANGULAR MOMENTUM FLUX IN 
ADM COORDINATES 



At this juncture let us remind where we are eventually headed. We are interested in the phasing 
of the binaries moving in quasi-eccentric orbits and in the first instance, as for quasi-circular orbits, 
we work in the adiabatic approximation. In this limit the radiation time scale would be much 
longer than the orbital time scale and consequently we would require an averaged description of 



6 The results can also be derived directly in the general radiation gauge by explicit use of the mass quadrupole 
expression in the general radiation gauge Eq. (13. 1U and the 2.5PN equation of motion terms in the general radiation 
gauge Eq. ( 13.7b . Alternatively, it can also be obtained by a direct transformation of the flux expressions in SH gauge 
to general gauge using the shift (13 .9b and the equations (6. 1)— (6.3) from I20I1 for the changes in r, r and v 2 , together 
with a similar relation for (p. We have verified that these various methods lead to the same results. 
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the radiation reaction over an orbital period; so we average the flux over one orbit. To deal with 
3PN radiation reaction (corresponding formally to 5.5PN terms in the equation of motion) we 
require the description of the motion to be 3PN accurate. Such a description is available from 



the work of Memmesheimer, Gopakumar and Schafer [12811 based on the known 3PN equations of 
motion [[[], Q, 0, @], and we shall use this generalized 3PN quasi-Keplerian (QK) representation 
of the motion in ADM coordinates to average our angular momentum flux. The details of the 
expressions we use are available in Ref. [280 and we have recast these into forms better suited for 
the present context in Ref. ll20ll (see Sec. VII there). 

Since the quasi-Keplerian orbit is planar, to quantify the evolution of the orbital elements under 
radiation reaction we only need to average the magnitude of the angular momentum flux over an 
orbit. The computation thus becomes a generalisation of our earlier computation of the average of 
the energy flux in rf20ll and requires similar intermediates. Using the QK representation of the orbit 
in ADM coordinates and the instantaneous angular momentum flux in ADM coordinates obtained 
in Sec. [TTTl one transforms the expression for the norm of the angular momentum flux @ inst (r, r, (p) 
into another expression @ iBSt (E, h, e t , u) depending on the 3PN conserved orbital energy E, on the 
norm of the 3PN conserved angular momentum J = \J'\ as rescaled by Gm (thus h = J/Gm), on 
that particular choice of eccentricity e t which parametrizes the (3PN-generalized) Kepler equation, 
and on the eccentric anomaly u; see Eqs. (7.1)-(7.4) in rf20ll . Like for the case of the energy flux 
we find that the general structure of the angular momentum flux in terms of these variables reads 
up to 3PN order, 



dw 1 



a N (E,h) sinw ln(l-e,cosw) 
+ f3 N (E, h) + y N (E, h) 



(1 - ^cosw)^ (1 - e, cos u) N (1 - e t cos u) 1 



(4.1) 



For later convenience we have factored out du/d£, where £ = n(t - t P ) denotes the mean anomaly, 
with the instant of passage to the periastron, and n = 2n/P is the mean motion, with P the 
orbital period. The coefficients a N , fi N and y N are explicit functions of the invariants E and h, 1 
and are straightforwardly deduced from the QK parametrization in the form of a PN series — see 
Eq. (4.12) of ll25ll for instance — but too long to be listed here. Rewriting the angular momentum 
flux using the generalized QK representation, the flux can be averaged over an orbit to order 3PN 



extending the results of Il25ll at 2PN. The orbital average is (setting f P = for definiteness) 



(^inst) 



1 r 1 r d£ 

- &mst(t)dt = — g imt (u)—du. (4.2) 
P Jo 2ttJ dw 



To compute it we use some integration formulas. First it is clear that the second term in (14.11) will 
not contribute, 



1 f 

2tt Jo 



sin u du 



(1 - ecos u) N 

Note that this term corresponds to the 2.5PN contribution in the angular momentum flux which is 



7 Notice that the structure of Eq. ( 14.1b requires and /3n to be functions of E and h. The angular momentum flux 
could be averaged without reducing it to this form, though in alternative forms the average of the equivalent of the 
7ai term may require another standard integral different from that given in Eq. ( 14.6b . 
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therefore seen not to contribute to the average. To get the first term we have the useful formula 
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(4.4) 



where the right-hand-side (RHS) is to be evaluated at the value y = 1 after the N- 1 differentiations. 
We have also an alternative formulation in terms of the Legendre polynomial Pn-i, 
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(4.5) 



Furthermore, we dispose of the more complicated formula Il20ll 
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(4.7) 



These formulas, and notably the third one (I4.6I) - (I4.7I) . permit one to display the result in a com- 
pletely closed form, and give 
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(^inst) = ^ 
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(4.8) 



The results can be expressed in terms of different choices of variables as for the energy flux 
in |[20ll . For reasons discussed there, we choose the pair consisting of (e t , x) with 



x = 



iGmu) v 



2/3 



(4.9) 



where the mean orbital frequency is co = (<p) = Kn, equal to the mean motion n times the peri- 
astron precession K = 1 + k (with k denoting the relativistic precession). The expression for the 
orbital averaged angular momentum flux can be finally written (in ADM coordinates) as 

<&nst> = ^c 2 mv 2 x 7/2 (k n + x <H 1PN + x 2 K 2m + x 3 -Kspn) , (4. 10) 



where the individual PN terms read as 
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This ADM-coordinates expression, similar to the energy flux case, does not contain the logarith- 
mic terms ln(r/r^), consistent with these being pure gauge dependent terms arising specifically in 
the SH coordinate system. However, we notice that there are still some other logarithmic terms 
involving the constant r , in the instantaneous part of the flux, even in ADM coordinates. Actually 
these terms are parametrized by x in ( |4.11d| ) where 



Xq 



( j i a 

c 2 r Q 



(4.12) 



Recall that r = c r is an arbitrary length scale introduced in the general MPM formalism (to 
regularize Ultra- Violet divergences), which then appears in the definition of the source multipole 
moments starting explicitly at the 3PN order. This is what leads to the In r dependence of the 
instantaneous terms in the angular momentum flux. It is known [26[ 55] that the lnr terms will 
be cancelled by similar terms when we add up all the non-linear multipole interactions (mainly 
tails) constituting the radiative multipole moments observable at infinity. This has been checked 
at the 3PN order for compact binaries in the circular orbit case [0], where the contribution due to 
tails-of-tails effectively removes these unphysical lnr terms. In Ref. ll20ll . this cancellation was 
proved in the case of the 3PN energy flux for general non-circular orbits. We shall also recover 
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this necessary cancellation below for the angular momentum flux for general orbits. 

As a check of the algebra, we take the circular orbit limit of the averaged instantaneous angular 
momentum flux, say (Qmst)o, as given by (I4.10I) - (I4.1 II) with e t = 0, and get 



<0mst>0 = yc 2 mv 2 x 7/2 



1247 35 \ / 44711 9271 65 , 
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94403 
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775 
324 1 



(4.13) 



This is in perfect agreement with the averaged instantaneous energy flux for circular orbits, say 
(^inst)o> m the sense that 

(4.14) 

where the proportionality factor is the binary's mean orbital frequency to = c 3 x 3/2 /Gm. For the 
reader's convenience we recall the 3PN energy flux (Timt) in Eqs. (16.61) — (16.7b below. 

We emphasize again that the expressions we provide for the orbital average of the angular mo- 
mentum flux and later the evolution of orbital elements are in the ADM coordinates; we can easily 
obtain the corresponding expressions in another coordinate system like the Modified Harmonic 



(MH) one H20L 12811 by transforming the eccentricity e, = ef DM used here to another eccentricity 
such as ef H . There is no need to transform the parameter x which is a gauge invariant. The relation 
linking ef™ and e™ u reads (see Eq. (8.21) in Bl) 
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(4.15) 



Since the transformation begins at 2PN order, the eccentricity in the brackets of the RHS of (14.151) 
can equivalently be replaced by e^ DM or e^ H . Furthermore only the instantaneous part needs to be 
transformed. The hereditary contributions retain the same form up to the 3PN order to which we 
are concerned at present. For the explicit form of the angular momentum flux in MH coordinates, 
see Appendix O 

Note also that the x parametrization that we use can easily be changed to the alternative 
parametrization by the variable 

Gmn 

f=— 3~. (4-16) 

which is often used in the literature and is also gauge invariant. In ADM coordinates we have 
(with e t = ef DM ) 
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(4.17) 
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V. HEREDITARY CONTRIBUTIONS TO THE 3PN ANGULAR MOMENTUM FLUX 



Having obtained all the instantaneous terms in the 3PN angular momentum flux (in averaged 
form), one must now turn attention to the hereditary contributions. As for the energy flux 112 ill 
it is not feasible to obtain closed-form results in the time domain for these contributions, and 
we adapt the strategy set out there by using a (discrete) Fourier decomposition. The details and 
notation are similar to those in 112 ill so we avoid repeating them and only quote the final results 
with a few intermediate steps. We however find a new aspect that arises in the case of the angular 
momentum flux with respect to the energy flux, namely the presence of a contribution due to the 
memory integral in the angular momentum flux for general systems. We shall prove that this 
memory contribution yields an interesting zero-frequency effect depending on the past history of 
the system (see Appendix lAl). 



A. Tail and tail-of-tail integrals at Newtonian order 

Most of the hereditary contributions will require only relative Newtonian precision. At Newto- 
nian order we can use the Fourier decomposition of the Newtonian (N) source multipole moments 
/[ N) and J^_\ given by 

+00 

— (p) 

p=-ca 
+00 

^ (p) 

p=—oo ' ' 

with discrete Fourier coefficients ( P )i^ N) and ( P )3' ( ^\ indexed by the integer p. Since the moments 
are real we have e.g. (~ P )I ( ^ ) = (p)^i where * is the complex conjugation. Here we denote the 
mean anomaly by I = n(t - t P ), with n = 2n/P being the mean motion and P the orbital period (we 
can choose the origin of time so that t P = 0). The latter Fourier decompositions are simple because 
the motion is periodic at Newtonian order since there is no relativistic precession, K = 1 + 0(c~ 2 ). 



We follow exactly the same method as in Ref. H21H . and express the tail, tail-of-tail and tail 
squared terms (12.9I )- (12.10| ), i.e. up to the 3PN order, and averaged over the mean anomaly £, in 
terms of the Fourier coefficients of the multipole moments. All these terms are Newtonian except 
the mass-type quadrupolar tail term given by the first term in (12.91 ) and which must be evaluated at 
1PN. For the Newtonian part of the mass quadrupole tail — a quadratic interaction oc G 2 between 
the mass quadrupole moment /y and the mass monopole M — we get 

( £taiK(N) = _ i 8* G^M_ Y(vn) 6 J (N) T (N) * (5 2) 

\»i Vass quad 1 <: „1() /S" U > f Ja f ka ' \ J -^> 

where the range of p's is set to correspond to positive frequencies only. We add a label (N) to make 
the distinction with the mass quadrupole tail we compute below at 1PN order. The remaining tail 
integrals involve the Newtonian mass octupole and current quadrupole moments and read 

A /~ , 2 A /T +cx> 

(^)mass oc, = -i | ^ S ijk |> nf J$ 1$ . (53a) 
0.5 C ^ {pY {p) 
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40 c " (?) (p) 

In Sec. IV El we shall provide the plots, computed numerically, for the relevant "enhancement" 
eccentricity-dependent factors associated with Eqs. (15.31) . since they do not admit a closed-form 
expression. 

In addition to the previous quadratic-order tails, we have also at the 3PN order the first cubic 
non-linear interactions between /y and two mass monopole factors M, namely the so-called tail-of- 
tail integral and the tail squared one, both being evaluated at Newtonian order. For the tail-of-tail 
part, averaged over an orbit, we get 

/<g tail(tailK _ • S <J M „ WnnY 7 7~ (N) T (N) * 

in 2 t ,2 57/ v 1246271 

x J _ _ 2 (ln(2/7 n r ) + c) + — (ln(2p n r ) + c) - ■ (5-4) 

Here C = 0.577 . . . denotes the Euler constant. Note that in contrast to the quadratic tails this 
expression involves some logarithms, and even some squares of logarithms. However some more 
logarithms are contained in the related contribution of tail squared, namely 

(&r'> = -i J~Tr *M {y + ^P"ro) * C - ji) 2 } , (5.5) 

and we can check that in fact the squares of logarithms cancel each other when we add together the 
two contributions (15.41) and (15.51) . Such cancellation is known to occur for general sources [26]. 
Hence we get for the sum 

z^taUctaiD+ctaiifx 8 G 3 M 2 ^ - (N) (N)t ( 2n 2 214/ x 1167611 

} > = - 5 — ^ 1[ P H) $ fa {— ~ 105 ( ln(2p H To) + C ) " ^9400" j • 

(5.6) 

This result still depends on the arbitrary time scale r . It will be important to trace out the fate 
of this constant and check that the complete angular momentum flux we obtain at the end is 
independent of To = r /c. 

It is worth recalling that the reduction of the above expressions to their simple form is made 
possible thanks to the following closed-form formulas which are applied to each of the Fourier 
components of the flux when performing the time average. For any cr denoting the product p n 
(which at the initial stage can be positive or negative), we have 

(5.7a) 



^sign((r) + i(ln(2|o-|T ) + C) 



I 



- <T R - V \ _ . e~ itTTR (n 2 
2r / cr 



dVe- itrV ln 2 \— — -I = i 



(5.7b) 



where sign(cr) and \cr\ denote the sign + and the absolute value of cr, and where C is the Euler 
constant. The rational fractions such as 11/12 present in the tail integrals are easily taken into 
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account by redefining the constant To in these formulas as To e~ nln for instance. For dealing with 
higher non-linear tails (occurring at 4PN order for instance), we would need further integration 
formulas involving higher powers of the logarithms. 



B. The mass quadrupole tail at 1PN order 

One hereditary contribution will require the relative 1PN order, namely the leading quadratic- 
order mass quadrupole tail. In this case the Fourier decomposition is more complicated and one 
must exploit the "doubly periodic" nature of the 1PN dynamics in the two variables I and A = K£, 
where K = 1 + k is the advance of the periastron per revolution. We can then write the following 
doubly periodic Fourier decomposition of the mass quadrupole moment at the 1PN order Il2in . 



j (p+m k) C 



(p,m) 



(5.8) 



The discrete Fourier coefficients ( p , m )Zij depend now on two integers: peZ and the "magnetic" 
number m e Z such that |m| < / = 2. We have ( _ p _ m) J,y = ( P , m )J* 7 - These Fourier coefficients are 
valid through 1PN order. Actually one can check that those for which m = +1 are zero in the case 
of the mass quadrupole moment at the 1PN order. 

We insert the Fourier decomposition (15.81) into the 1PN mass quadrupole tail [first term in 
Eq. (12.91 )1 and perform the average. The result in terms of the doubly-periodic Fourier coefficients 
at 1PN order reads 



mass quad 



4 G 2 M 



tail 



-1 - 



£ijk / n 7 (p + mk) 2 (p' + mk) 4 \p' - p + (m - m)k\ J j a I 



ka 



5 C , J (p,m) (p',m') 

p,p';m,m 

x (j(p+p'+(™™')W) ^ d y e -i(p'+^)»(^-v)j ln /Z|zZj + lij . (5 .9) 



The summations range from -oo to +oo for p and p', and from -2 to 2 for m and m'. The last 
two factors in (15.91) , i.e. the average over £ of an elementary complex exponential, and the Fourier 
transform of the tail integral, are both evaluated following 112 ill at first order in the relativistic 
advance of the periastron k which is a small 1PN quantity, k = 0{c~ 2 ). The ^-average factor reads 



/ A(j>+mk){\ _ f ^_ i(p+mk)C 



— k 



if p * 



(5.10) 



1 +inmk ifp = 



where we neglect some 2PN terms <9(c~ 4 ), and we have used the fact that m k «; 1 since we are 
in the limit where k — » (hence p + mk is never an integer unless k = 0). This result depends on 
whether p is zero or not, and is true for any integer m, except that when m = it becomes exact as 
there is no remainder term 0(c~ 4 ) in this case. The tail integral expanded to first order in k [i.e. up 
to some remainder <9(c 4 )] reads 
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(5.11) 

p l n 
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For the remaining integral in the RHS of (15.1 II) we apply the formula (I5.7al) . 

In this paper we do not give a more explicit form for the Fourier decomposition (15.91 ), which 
together with the two results (15.101) and (15.1 II) is well-defined. Indeed (15.91) is given by a very 
complicated expression which can only be handled using an algebraic computer program. Still it 
will remain to insert in that expression the explicit 1PN results for the mass quadrupole moment 
and the total mass M for eccentric binary orbits, and to re-express the series in terms of some 
elementary eccentricity-dependent "enhancement" functions which we shall evaluate numerically. 

C. The non-linear memory integral 

The memory contribution is defined, from Eq. (12.1 II ), by 

C emory = y 5 % £ dVJg(V) I2(V) . (5.12) 

In principle we should put brackets to indicate a STF projection acting on the indices k and a 
in the RHS, however the expression is manifestly symmetric with respect to those indices, and 
automatically trace-free thanks to the presence of e^. The orbital average reads 

<£~ y > = Y5 ?o I ^ J ^ dVlg\V) I%(V) . (5.13) 
We invert the two integration signs to re-write the latter expression as 

The contribution extending over the entire past (i.e. over -oo < V < 0) involves the orbital average 
of the third time derivative of the moment, i.e. 

(O-Jf ^(r*). (5 - 15) 

Since the quadrupole moment is Newtonian at this level of accuracy, the motion is periodic, and 
the quadrupole averages to zero. However one is not allowed to replace (TV ) = into the first 
term in the RHS of (15.141) because the argument neglects the evolution in the remote past of the 
Keplerian orbital elements by radiation reaction. In Appendix [A] we shall study the dependence 
over the binary's past history and find that it gives a contribution on the current dynamics in the 
form of a zero-frequency or DC effect. 

On the other hand the memory due to the recent past of the source [second term in the RHS 
of (15.141) 1 is easily seen to be zero on average. Indeed we evaluate the integral in the same way 
as was done to compute the orbital average of the instantaneous contributions in Sec. [IV] and find 
that it is made of a sum of elementary integrals only of the type (14.31) that are zero. Therefore we 
conclude that the memory contribution to the averaged angular momentum flux reduces to the DC 
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term due to the influence of the remote past of the source, 

<£™ y > = <^ C >- (5-16) 

The DC term has the structure of a Newtonian term and is obtained from a model of past orbital 
evolution for eccentric orbits in Eq. (IA12I) of Appendix lAl 8 

However, in this paper we have mostly in view the comparison with numerical simulations such 



as those in H45H . Numerical simulations start from initial conditions which are for the moment 
always set at some very recent instant. This means that the initial eccentricity e x is comparable 
to the current one e§ (using the notation of Appendix [A]). In this case one can neglect the DC 
contribution so that 

(g~ vy ) = 0. (5.17) 

In the present paper we adopt the result (15.171) appropriate for short-lived binary systems, but 
keep in mind the possible influence from the past of the non-linear memory DC term computed in 
Appendix lAl 

D. Definition of the eccentricity enhancement factors 

We shall now closely follow the numerical calculation of the hereditary terms in the energy 



flux [21]. We shall present here only the definitions we use and directly the results we obtain 



from those definitions. We refer to 112 lh for more details. The source multipole moments (in the 



center-of-mass frame) at Newtonian order read 

If = tisiv) 

Jf\ = ns l {v)x <L ' 2 s i '- l>ah x a v\ (5.18b) 



lf ] = ns,(v)x <L> , (5.18a) 



and involve the following function of the symmetric mass ratio v = /i/m, 9 

s l (v)=X 1 z 1 +(-) l X[- 1 , (5.19) 

where X, = and X 2 = Next, we rescale the source moments in 

an appropriate way and introduce the dimensionless moments I L and J h _\ by 

/[ N) = iua l s,(v)I L , (5.20a) 
J ( ^\ = fta'n s,(v) J L -\ , (5 .20b) 

where a is the semi-major axis and n = 2n/P is the mean motion (such that Kepler's law n 2 a 3 = 
Gm holds at Newtonian order). 

We now define a set of "enhancement" functions of the eccentricity e of the orbit (at Newtonian 
order) by means of the Fourier components of the rescaled moments. Such enhancement functions 
will exactly parallel similar functions valid in the case of the energy flux [21], and we shall adopt 



for these exactly the same notation as in 112 111 except that we add a tilde on these functions to 



8 Recently the DC non-linear memory terms in the gravitational-wave polarizations have been computed to post- 
Newtonian order in the case of quasi-circular binary orbits Jfi? 



9 Alternative forms for this function can be found in Refs. 15 ll 15211 
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distinguish them from the functions parametrizing the energy flux. We thus pose 



-t-oo 

/« - -r* L >I>» 5 IP 1™' ■ <5 - 21a) 

p=l 

+00 

He) = ~ e ijk U J] P 6 I fa II , (5.21b) 

lo (p) (p) 

P =i 

9ft ' + °° 

*° - ~ 16453 <5 ' 21C) 

yie) = -HeijkLtYyfjaJZ, (5.21d) 
U cp) (p) 

+00 

= ~ U J] p 1 h & > ( 5 - 21e ) 

j2 ^—f (p) (p) 

p=l 

. +00 

^ 2 ' 7ta <5 - 21f) 



Like for the definitions adopted in 112 111 all the latter tilde functions but one are chosen in such 
a way that they tend to one in the circular orbit limit, when e — > 0. The notable exception is 
jjf(e) which vanishes in this limit. This is easily checked since in the circular orbit limit (and at 
Newtonian order) the quadrupole moment possesses only the harmonic for which p = 2, and 
consequently the log-term in^(e) — Eq. (I5.21fl) — becomes zero. Most of these functions will 
not admit any algebraic closed-form expression, and we shall leave them in the form of Fourier 
series to be evaluated numerically. However, as we shall now see, two functions can be computed 
algebraically, namely f(e) and F(e). 10 

The first function parametrizes the Newtonian part of the averaged angular momentum flux, 

<£ (N) > = fc 2 vW' 2 />), (5.22) 

where x is defined by (14.91) and reduces at Newtonian order to G m/(ac 2 ). Thus, f(e) is the Peters 
& Mathews [22, 29)] function, which admits an algebraically closed-form expression which is used 



in the timing of the binary pulsar PSR 1913+16 16511 . and given by 

1 + 7 -e 2 

to- o3v ■ ( " 3) 

in agreement with the Newtonian part of our earlier result (14.10l) - (14.1 II) . 

On the other hand the function F(e) is the analogue of the function F(e) in rf2~0l and will partly 
parametrize the tails-of-tails, where it will have in factor a contribution depending on the arbitrary 
constant scale r . That dependence on r is the same as in the instantaneous part of the flux, as 
given by Eqs. (14.10l) - (14.1 II) . Such a specific dependence of the hereditary terms on r will just be 
appropriate to exactly cancel out the In r in the total angular momentum flux. The function F(e) 



We reserve Latin names for algebraic closed-form functions, and Greek names for numerically generated ones. 



25 



is given by 
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1 + 229 2 , 327 4 , 69. 6 
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(1 - e 2 ) 5 

The analogous function F(e) in the energy flux is recalled in Eq. (16.91) below. 

From the previous definitions we can express the quadratic tails at Newtonian order as 



(5.24) 
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(5.25a) 
(5.25b) 
(5.25c) 



In the Newtonian mass-quadrupole tail (I5.25al) we recognize in particular the coefficient 4-n com- 
puted analytically in Ref. [|24ll (recall that i£(0) = 1). The function <p(e) has already been computed 
numerically from its Fourier series (I5.21bl) in Ref. Il24n . Notice that these formulas and similar 
formulas below take exactly the same form (i.e. with the same coefficients) as the corresponding 
formulas valid in the case of the energy flux [HU]. The reason of course is that in the circular-orbit 
limit e —* the energy and angular momentum fluxes are proportional, and related by Eq. (I4.14I ). 
However, for non-zero eccentricities, the enhancement functions will differ from the correspond- 
ing functions in the energy flux, and this is why we add a tilde on them. 

In a similar way, with the above definitions we get the sum of contributions from tails-of-tails 
and squared-tails as 
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X(e)\. (5.26) 



The circular-orbit limit of this result is immediately read off and seen to agree with the previous 
result in Refs. fl,H. 

Finally we provide the mass quadrupole tail at relative 1PN order. Its computation is much more 
involved because the Fourier series (15.91) contains several summations, and depends on intermedi- 
ate results (15.91) and (15.1 II) . The computation is based on the known 1PN relativistic corrections 
in the mass quadrupole moment and the total mass M, which are given in Eqs. (5.16)— (5.17) 
in 112 ill . The result is of the type 



<£? tml >mass quad = y c* Vmx^\ 4n ip{e t ) + nx 
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(5.27) 



which defines two new enhancement functions a and 9 which are the analogues of the functions 
a and 9 in the energy flux 112 ill . Since a and 6 are given by some very complicated Fourier series 
(handled with Mathematica) — instead of relatively simple ones like in (15.211) — we shall directly 
compute them numerically using the same method as in ||2~Hl. Notice that since we are at the 1PN 
order we must be specific about which definition of eccentricity we use; here we adopt the time 
eccentricity denoted e, which enters the generalized Kepler equation at 1PN order 0211 . On the 
other hand, the 1PN corrections arising from the parameter x [see (14.91 )1 are evidently crucial in 
the result (they include the 1PN periastron advance k). 
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For the final presentation it is convenient to redefine in a minor way our elementary enhance- 
ment functions. Let us choose (still paralleling Ref. 112 IIP 
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(5.28a) 
(5.28b) 
(5.28c) 



Considering thus the 1.5PN and 2.5PN terms, composed of tails, and the 3PN terms, composed of 
the tails-of-tails and the squared-tails, we get the final form of the total hereditary contribution to 
the averaged angular momentum flux (12.81) : (£^ ered ) = L, (^hered) with 
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(5.29) 



For circular orbits this angular momentum flux is in agreement with the energy flux ('Fhered)© 
computed in [20], since we have <!Fhered)o = co (^hered)o- 



E. Computation of the enhancement functions 

The eccentricity-dependent enhancement functions we use in (15.291 ), computed numerically, 
are now provided in the form of plots. The numerical computation has been described in the case 



of the energy flux 112 111 and we adapt the same for the angular momentum flux. 

Essentially, the fitting procedure to obtain the Fourier coefficients of the Newtonian source 
moments (or their periodic part in the mean anomaly I for 1PN accurate moments) in terms of I 
can be implemented either starting with the basic multipole moments themselves or the leading 
time derivatives appearing in each of the terms in the angular momentum flux. The latter method is 
known to improve the numerical convergence of the final sum because one deals with lower order 
time-derivatives of the basic functions. However, here we have followed the former method which 
has the advantage that the fitting function is much more simple for the basic source moments than 
for their multi-time-derivatives and thus takes less time and is also less prone to errors. Proceeding 
in this way also provides another check on the energy flux calculation as we have reproduced 
the results of Ref. 1I21I1 using this alternative choice. At the Newtonian order it is in fact more 
efficient to make use of the well-known Fourier decomposition of the Keplerian motion to compute 
the Fourier coefficients. Using this we can derive the components of the multipole moments (at 
Newtonian order) as series of combinations of Bessel functions. It is then simple to compute 
numerically the associated Newtonian enhancement functions [namely the functions <p(e), /3(e), 
y(e) and^e)]. On the other hand, for the Newtonian tail terms, we could proceed exactly in 
the same way as for the 1PN term, following the various steps and evaluating numerically the 
functions. We have verified that both methods agree well. The above procedure is quite general, 
and provides a method which could be extended to higher PN orders. 

We present in Figs. [T|-[3] the functions which permit to define the hereditary part of the angular 
momentum flux (I5.29I ). We recall that these functions are such that they reduce to one in the 



27 




FIG. 1: Enhancement function Cp(e) in the angular momentum flux at 1.5PN order. 
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FIG. 2: Enhancement functions \fr(e) and 1(e) in the angular momentum flux at 2.5PN order. 



circular-orbit limit e — » 0. To facilitate the comparison with the results of numerical relativity [45] 
or use in data- analysis applications we provide also some numerical tables for these functions in 
Appendix El 



VI. EVOLUTION OF ORBITAL ELEMENTS UNDER 3PN RADIATION REACTION 



The most important application of the 3PN angular momentum flux is to calculate, using also 
the energy flux Il20n . how the binary's orbital elements evolve under 3PN gravitational radiation 
reaction. We shall compute the time evolution of the mean motion n, the periastron precession k, 
the mean orbital frequency co = nK, the semi-major axis a r , and the time eccentricity e t . 



A. General method 



We start with the 3PN expressions for n, k,co, a r and e t in terms of the conserved energy E and 
conserved angular momentum J of the orbit um- In an attempt to simplify some expressions we 
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often employ, in place of the energy E and angular momentum J, the dimensionless variables I120ll 

(6.1a) 
(6.1b) 



s = 



2E 
2EJ 1 



J = ~ 



(Gin) 2 



Since e = 0(c~ 2 ) this energy variable can be viewed as a book-keeping parameter labeling the 
successive PN orders. Recall that the semi-major axis a, and the eccentricity e, depend on the 
coordinate system, but that the mean motion n and periastron precession k are gauge invariant, i.e. 
take the same expressions in ADM and, say, modified harmonic coordinates. Of course this is also 
true of the mean orbital frequency co and hence of the parameter x. The expressions we give for e t 
and a r are valid in ADM coordinates. We have 
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The procedure to compute the evolution of the orbital elements under gravitational radiation- 
reaction is straightforward but lengthy. Differentiating the orbital elements with respect to time, 
and using the heuristic balance equations, we equate the decreases of energy and angular momen- 
tum to the corresponding averaged fluxes, and obtain the (secular) rate of change of the orbital el- 
ements. This extends earlier analyses at previous PN orders: Newtonian Q22J] , 1PN order Il23l.l30ll. 
1.5PN order [24], 31] and 2PN [25, 37]. Taking the example of the mean motion we have 



dn dn dE dn dJ 
~dt ~ lm~dt + ~dJ~dt 



(6.3) 



The usual (heuristically derived) balance equations for energy fi(dE/dt) = —(T) and angular 
momentum // (dJ/dt) = —{§), where the fluxes are known up to 3PN order, give the 3PN evolution 
equation averaged over one orbit, 



dt fi 
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For eccentric orbits we recall that this gives only the slow secular evolution under gravitational 
radiation reaction. The complete evolution includes superimposed on this a fast but smaller peri- 
odic oscillation on the orbital time scale which can be conveniently computed using a two-scale 
decomposition following IU9L 13711 . 

To express the final result of (16.41 ) in terms of x and e t , the variables in terms of which we have 
represented the energy and angular momentum fluxes, we will require the expressions for e and j 
in terms of x and e t . We have, 
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B. Recalls for the energy flux 



To proceed further we also need the energy flux and so we recapitulate here the results from ||2 
2 111 . The instantaneous part of the energy flux reads 

32 c 5 / \ 
< ^inst) = y — v 2 x 5 \I N + xIi PN + x 2 J 2PN + x 3 J 3 p N j , (6.6) 
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where the coefficients (in ADM coordinates) are 
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The hereditary part of the energy flux is given by 
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With the exception of F that can be given in a closed analytic form, 
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the untilded enhancement functions tp,ifr, ... (differing for non-zero eccentricities from the tilded 
ones <p, iff, . . . in the angular momentum flux) are computed numerically in bill ; in the Appendix|B] 
below we provide the numerical tables of these functions. 



C. Instantaneous contributions 



It is clear that since the evolution of the orbital elements is linear in {T) and {Q) one can 
separate out the contributions due to the instantaneous and hereditary components in the fluxes. 

Starting with the instantaneous terms, we find for the evolution of the orbital elements 
{n,k,co,a n e t } n a PN structure with coefficients depending on E and J or, alternatively, on the 
frequency-related parameter x and the time eccentricity e t . Choosing an eccentricity parameter 
like e t as one of the basic independent variables has the advantage of yielding more usual-looking 
formulas, e.g. by recovering at the lowest order the Peters-Mathews enhancement function. How- 
ever it has the disadvantage of depending on the employed gauge. For some purposes it is better 
to use a pair of gauge invariant variables like (E, J) or its variant (e, j). Another interesting choice 
for apair of gauge-independent variables is (jc, k) or alternatively (x, i) with i = 3x/k as was used 
in Il20n . Here we present the results in terms of the pair of parameters (x, e t ); if necessary it is 
straightforward to transform the results into a different set of parameters like (x, k) or (e, j). 

We start with the time evolution of the mean motion n = In I P. Since n depends on the angular 
momentum J only from the 2PN level [see (I6.2al) l, the angular momentum flux (Q) will only be 
needed in this case up to 1PN order (while the energy flux is needed with full 3PN accuracy). We 



Of course this is a redundant set of orbital elements; for instance co = n (1 + K). 
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present the result in the form 
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Next, the evolution of the periastron precession is 
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Note that there is no dependence on the scale x in this expression; the reason is that the 3PN coef- 
ficient is actually 2PN relatively to the dominant term. From these above results we immediately 
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deduce the time evolution of the mean orbital frequency to = Kn as 
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For the semi-major axis a r , again we need the 3PN energy flux but only the 1PN angular momen- 
tum flux, and get 
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The evolution of the eccentricity e, is the only one to require both the energy flux {T) and angular 
momentum flux (Q) with full 3PN accuracy. We obtain 
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The leading order Newtonian term is in agreement with the work of Peters I122I1 . The match to 
earlier results for the evolution of orbital elements includes also 1PN lf3Qh and 2PN |25] orders. 
The explicit expressions for the evolution of the orbital elements in modified harmonic coordinates 
are provided in Appendix O 



D. Hereditary contributions 

The hereditary contribution to the flux begins at relative 1 .5PN order and consequently the 1PN 
quasi-Keplerian representation, given by the truncation of Eqs. (16.21) at the 1PN order, suffices for 



39 



this analysis. Namely, 
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We also need the 1PN accurate expressions for s and j, which we express in terms of e t and x: 
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The results for the hereditary parts (!Fh er ed) and (inhered) depended on some numerically- 
computed enhancements functions of the eccentricity: ip, if/, g, ... in the energy flux, and Sp, 
iff, . . . in the angular momentum flux. Obviously the evolution of orbital elements will involve 
some linear combinations of cp, if/, . . . and their tilde analogues [see e.g. (16.41) 1. We thus have 
to introduce some new enhancement functions to parametrize the time evolutions of the various 
orbital elements {n, k, a>, a n e t ). We pose 12 
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To avoid heavy notation and since we deal with only the time eccentricity e t , the corresponding enhancement 
functions <p e , if/ e , g e , K e and F e are labelled by e rather than e,. 
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Note that £ n (e) = g a (e) = £y(e). We have also a function known analytically, 
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I these new functions ol eccentricity reduce to one in the circular orbit limit e — » 0. 
We now list the relevant 1PN accurate expressions for the hereditary part of the evolution of 
the orbital elements, as computed from the hereditary energy and angular momentum fluxes ( 16.8I) 
and (I5.29I) . However we notice that, for what concerns the hereditary part, the evolutions of n and 
a,- depend only of the energy flux (I6.8I) up to the accuracy we need (which is 1 .5PN relative order); 
similarly the evolution of k depends only on the angular momentum flux. We find 
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The hereditary part of the evolution of e t depends on both the energy and angular momentum 
fluxes at 1 .5PN order, and we find 
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FIG. 4: Enhancement functions ip(e) and (p e (e) in the evolutions of n and e t at 1.5PN order. 




FIG. 5: Enhancement functions if/ n (e) and i^ n {e) in the evolution of n at 2.5PN order. Note the similarity of 
tfr n (e) with ifr(e) in Fig. [2] 
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All the numerical functions of eccentricity are provided in Figs. |4]-[T01 Numerical tables are given 
in Appendix |B] (see Tables I-IV). 

VII. LIMITING FORMS FOR SMALL ECCENTRICITY 



In some cases one may have prior information about the smallness of the eccentricity and in this 
case one may only need the leading corrections when the eccentricity parameter e, — » 0. The main 
problem we face is to treat the hereditary parts because we could compute them only numerically 
for general eccentricities. However it was shown in 112 ill how to obtain analytically the leading 
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CD 




FIG. 7: Enhancement functions if/ e (e) and £ e {e) in the evolution of e t at 2.5PN order. Note the similarity of 
i// e (e) with if/{e) in Fig. [2] 

corrections [neglecting 0{e A t )] for the enhancement factors in the case of the energy flux. However, 
here we need to push the accuracy of these results to the next order [neglecting 0(ef )] in order to 
evaluate the leading-order correction in the enhancements functions present in the evolution of the 
orbital eccentricity, Eq. (16.251) : this can be checked from the explicit expressions of these functions 
as given in (16.221) . We find 
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FIG. 8: Enhancement functions K e (e) and F e (e) in the evolution of e t at 3PN order. 




FIG. 9: Enhancement functions <pk(e) and iff^e) in the evolution of k and co at 1.5PN and 2.5PN orders. 
Note the similarity of \j/ w (e) with if/(e) in Fig. [2j 
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Similar results are obtained in the case of the angular momentum flux, i.e. for the tilded enhance- 
ment factors, which read 
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FIG. 10: Enhancement function if/ a (e) in the evolution of a r at 2.5PN order. The other function £ a (e) at 
2.5PN order is the same as C, n {e) and hence is not plotted. Similarly at 3PN order it involves the same 
functions /c(e) and F(e) as for n. Note the similarity of with \J/{e) in Fig. [2} 
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Using (17.11) and (17.21) we then get the leading-order corrections in the functions parametrizing the 
evolution of orbital elements and defined by Eqs. (I6.22M6.23I ), 

94115 

Me) = 1 + Yfs9§ e +0(e 4 ) , (7.3a) 

Ue) = l + ^e 2 +0(e 4 ), (7.3b) 
257 

WtGO = l + —e 2 +0(e 4 ) , (7.3c) 

Me) = l- 5 -^-e 2 + 0(e 4 ), (7.3d) 

= 1 + ^T e2 + ( e4 )' (73e) 

9 1 79Q 

= l + ^^^ + O^), (7.3f) 



45 



K e (e) = 1 + 



1 637 339 503 + 2 135 608 720 In 2 - 663 415 515 In 3 



j ^ , 14023 2 
F e (e) = 1 + — — e 2 + 
1538 

82775 , 
Me) = 1 - e 1 + 



4159 



179 578 418-97 909 280 In 2 + 98 283 780 In 3 

o(*) , 



] je 2 +0(e 4 ), (7.3i) 

(7.3j) 
(7.3k) 
(7.31) 



To finish, we provide the complete results (composed of both instantaneous and hereditary con- 
tributions) valid up to first order in e 2 . The total fluxes of energy and angular momentum read 13 
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13 Consistent with the accuracy which was needed in Eqs. ( I7.U -( |7T2] |. we recall that the evolution of the orbital 
eccentricity to the level ej [as given by ( I7.6eb below] requires these expressions to be accurate up to e\ . 
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We have compared the e 2 terms of the angular momentum flux expression above with the result of 
black-hole perturbation theory given in Eq. (181) of Ref. [66J]. Similar to the case of the energy 
flux [20], we find that our result in the test-mass limit v — » matches with the perturbation result 
with a transformation of eccentricity given by: 
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Here the eccentricity e t is the one appearing in (I7.4bl) . i.e. in ADM coordinates, while e is the 
eccentricity used in perturbation theory. The relation (17.51) is the same as the one we found for the 
case of the energy flux, as it must surely be: indeed see Eq. (9.6) in rf20ll where we must take into 
account the change between the modified harmonic (MH) coordinates used there and the ADM 
coordinates used here. 14 

The evolution of orbital elements is 
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From Eq. ( 14.151 ) we have ef H = e^ DM [l - |x 2 - 4x 3 + 0(v)] in the test-mass limit for small eccentricity O(ef). 
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(7.6e) 



We find that in the limit when e, —* 0, the evolution of the orbital frequency, i.e. (dco/dt) Q , 
reproduces the known result for circular orbits at 3PN order as given in Q8L |9Q . 
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APPENDIX A: THE NON-LINEAR MEMORY DC TERM 

Our aim is to compute the zero-frequency part (or DC part) of the non-linear memory integral 
found in Sec. IV CI namely 

C em ° r Vo) = | ^ s ijk I$(t ) £ dtlg\t) />) . (Al) 

For convenience in this Appendix we denote the current time of the observer by t Q (formerly 
denoted T R ) and the earlier time over which the memory is integrated by t (formerly V). In Eq. (IA1I) 
we shall suppose that the system was formed at some initial instant t\ on some very eccentric orbit 
with initial eccentricity e\ close to 1, but strictly less than 1. Then the system will evolve by 
radiation reaction until reaching the current eccentricity eo such that e Q e\ < 1. 

At this order of approximation we can replace in (IA1I) the quadrupole moment and its time 
derivatives by Newtonian values. We express the result in terms of the current orbital separation 
r = r(to), radial velocity r and instantaneous orbital frequency <p - The result will be an integral 
extending on the corresponding variables r, r and <p for the orbit at any instant in the past [recall 
our notation after Eq. (13.11) ; here e.g. r = r(t)\ For convenience we set the origin of the orbital 
phase (f (or true anomaly) at the current binary's separation, i.e. we pose cp = 0. A straightforward 
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calculation gives the norm of the memory integral as 



^memory 64 G 5 mV f° J 17 f 2 r (f A „ / (ft 1 r Q f \ . „ 

^memory = ^^^ df ^ I 2 _ 4 ^ _ sm2 ^ 

105 c w r J h \\ r 4 r r 3 / \ r 3 4r r 4 / 

We then replace r, f and <p in this integral by the explicit solution for Keplerian motion, 

a(l - e 2 ) 

r= ~i ' 

1 + e cos (p 



(A2) 



(A3a) 



r= A /— : ^-esin^, (A3b) 

a(l - e l ) 

-^—^(1+ecos^) . (A3c) 

The resulting integral is composed of many terms corresponding to different harmonics of the 
Keplerian motion. Thus, we find that the structure of the memory term is of the type 



^memory = £ P df/^e)^, (A4) 



where the /„'s depend on the semi-major axis a{t) and eccentricity e{t) of the orbit in the past (and 
depend also on the current values r , r , (p ). Here (p(t) is the orbital phase which is oscillating at 
the orbital period. In contrast, a(t) and e(t) remain approximately constant during one period, but 
slowly evolve by radiation reaction during all the binary's past evolution. Now we expect that the 
oscillations of the orbital phase, due to the sequence of many orbital cycles in the entire life of the 
system, will more or less cancel each other in the memory integral, so that the true post-Newtonian 
order of the oscillating terms will be simply given by the power of 1/c they carry. The oscillating 
terms, having n + 0, are thus expected to have their normal 2.5PN order, 15 and we have shown in 
Sec. IV CI that these terms do not contribute to the averaged angular momentum flux. 

With this, now there remains the purely zero-frequency or DC component of the memory, 
corresponding to n = 0, which is given by 

& DC = f°dtMa,e). (A5) 

This DC integral is cumulative because it extends over some steadily increasing kernel and a priori 
exhibits a strong dependence on the past (i.e. when t\ — » -oo). The explicit calculation gives 



16 G 6 mV y p e 2 

c 10 r J tl a\\-e 2 f 



, 20 2 19 4 

1 + — e + e 

21 336 



(A6) 



We now evaluate the DC term by inserting the secular evolution of orbital elements to Newtonian 
order consistently. We first convert the time integration into an integration over eccentricity by 



15 This expectation could be justified in more details using calculations similar to the ones in Sec. 4 of 14911 



50 



using the Peters Il22ll formula 



304 G 3 m 3 ve / 121 



* = -~ij wi-sM 1+ ^ ] - (A7) 

This is the dominant term in Eqs. (|6.18I )-( I6.19I ); for simplicity we ignore the averaging procedure. 
Calling e\ the initial eccentricity at the instant of formation of the binary system (e\ > e ), we 
obtain from Eq. (IA6I) 

3_GWn f , e 1 + 2 £e 2 + ^ 
* ~ 19 r J eo a(l-e 2 ) 5 / 2 1 + %e 2 " 1 j 

We discover here the memory effect: indeed, the real post-Newtonian order of the DC term (IA8D 
is found to be 1 /c 5 instead of the formal order 1 /c 10 exhibited in Eqs. (IA1I) - (IA2I) . which means 
that it occurs at Newtonian level in the angular momentum flux instead of the 2.5PN order. This 
increase by a factor 0(c 5 ), which corresponds to the inverse of the dominant order of radiation 
reaction, is clearly due to the cumulative integration over the past of the zero-frequency mode. 

We can simplify the result by using the relation for a{e) deduced from the evolution equations 
for e(f) and a(t). Although the exact expression is known, we rather employ for simplicity an 



approximate relation given in Ref. II22I1 . namely 

a 1 - e- Q i e 



I -el ' - ^ 12/19 



a 1 - e 2 \eo 



(A9) 



This finally yields 



c dc _ 3 Ij m v <po e / e 1 + 2 i g + 336 g rAim 

19 r flo(l - e 2 Q ) J eo ^ (1 - e 2 f> 2 1 + ^e 2 ' K ] 

We can see the importance of the dependence over the past by the fact that the integral diverges 
when the initial eccentricity e\ of the orbit approaches one [this would also be true had we used 
a more exact relation for a(e)]. For a binary system born on a very elongated elliptic orbit close 
to a parabolic one, the DC term is dominated by the initial evolution when e\ —* 1 and reads 
approximately 

--— — —-—7 — t,^=- (AH) 

119 c 5 r a (l-e 2 ) 



Finally we are interested in this paper in the averaged angular momentum flux. The formula (1A10I) . 
averaged over the current orbit, gives immediately 

/^>DC\ 3 2 2 7/2 e J j e 21 c 336 c /A11\ 

® > = - T 9 myC *° (T^fl ^ O-e^ l + ^ ■ (A12) 

This indeed appears, comparing with (I4.10l) - (l4.1 II) . as a Newtonian-like term. This term will mod- 
ify the enhancement factor at Newtonian order, namely (1 + |e^)/(l — e 2 ,) 2 , by an extra contribution 
coming from the binary's earlier evolution. However we expect that this could be important only 
for very long-lived binary systems having started on a nearly parabolic orbit. For such systems the 
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memory-induced modification of the Newtonian enhancement factor could be used to obtain an 
improved model of evolution of the orbital elements. 



APPENDIX B: TABLES OF NUMERICAL RESULTS 



To facilitate the quantitative comparison of the PN results with high precision numerical com- 
putations of the inspiral and merger of eccentric binaries |45[|, we provide the numerical tables of 
all relevant enhancement functions in this Appendix. 
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TABLE I: Tables for the numerical enhancement functions appearing in the expression of the 3PN hereditary 
energy flux in Eq. (16.81) . 
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TABLE II: Tables for the numerical enhancement functions appearing in the expression of the 3PN heredi- 
tary angular momentum flux in Eq. ( 15.291 ). 



APPENDIX C: ANGULAR MOMENTUM FLUX AND EVOLUTION OF ORBITAL ELE- 
MENTS IN MODIFIED HARMONIC COORDINATES. 

In this paper we have first obtained the angular momentum flux in the standard harmonic coor- 
dinates [Eq. (13.41 )1 - in terms of r,r,tp- and transformed it to expressions in the ADM coordinates 
rEq. (l3.5l) l. All subsequent formulas of the averaged flux and evolution of the orbital elements 
under gravitational radiation reaction refer only to ADM coordinates. Recent studies of binaries 
moving in elliptical orbits also employ the modified harmonic coordinates and for the convenience 
of such investigations we provide in this appendix explicit forms of the important equations in 
modified harmonic coordinates. 

Following the prescription outlined in Sec. VI A of Ref. ll20ll (see Eqs. (6.2) and (6.3) there), 
one can compute the angular momentum flux in modified harmonic (MH) coordinates starting 
from the corresponding expression in the standard harmonic (SH) coordinates in Eq. (13.41) . The 
difference between the standard harmonic and modified harmonic coordinate is a 3PN term. The 
final result in MH coordinates may be written as Q M h = Q S h + S^Q, where 

704 G 5 mV \r 2 ( v 2 3r 2 Gm\ / r Y| 
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TABLE III: Tables for the numerical enhancement functions appearing in the 3PN hereditary part of the 
evolution of orbital elements n, k and a> in Eqs. ( 16.241 ). Recall that = 



As expected, in the expression for the angular momentum flux in MH coordinates the lnfr/r^j 
terms will be gauged away. 

We next provide the expressions for the angular momentum flux averaged over an orbit and the 
secular evolution of various orbital elements in the modified harmonic coordinates. Notice that, in 
the 3PN accurate expressions, only the instantaneous terms at 2PN and 3PN will be different from 
the corresponding ADM expressions given earlier since the difference between the ADM and MH 
coordinates are at order 2PN and higher. Thus all the hereditary terms (which start at 1 .5PN) will 
be the same in the ADM and MH coordinates up to 3PN. Starting from the ADM expressions, 
one can obtain the corresponding results for various quantities by the sole transformation of the 
eccentricity parameter e t as given in Eq. (14.151) . Listed below are the corresponding expressions at 
2PN and 3PN orders in which e, now stands for ef 111 . As is obvious, in the equation for evolution 
of K, only the 3PN term will be different. 

^, MH If 135431 11287 260 , / 598435 9497 1546 ,\ , 

« ■ ^eV\-— + — Y+ -^ + V— + l^ V+ —T> 
30271 106381 „\ „ / 30505 2201 1519 
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TABLE IV: Tables for the numerical enhancement functions appearing in the 3PN hereditary part of the 
evolution of orbital elements a r and e t in Eqs. (|6.24| )-( [6.25l ). Recall that g a (e) = £ n (e). 
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We conclude this appendix with two useful relations, the analogues of Eqs. (14.151) and (14.171) 
but expressed in MH coordinates. We have, 
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